############
#	Lab13
############



# set up
compet = function(N1, N2, K1 = 150, K2 = 100, a12 = 1, a21 = 1){
	#K1, K2: carrying capacities
	#a12, a21: competition coefficients
	#n1post is the post-interaction Sp1 density
	n1post = N1*(K1 - N1 - a12*N2)/K1
	#n2post is the post-interaction Sp2 density
	n2post = N2*(K2 - N2 - a21*N1)/K2
	#the two vectors of results are stored in a 'list'
	res = list(N1 = n1post, N2 = n2post)
	#the list is passed out of this function
	return(res)
} #end function


#1. define T 
T= 100
#2. two population vectors
N1 = N2 = rep(NA, T)
#3. Set initial population values
N1[1] = N2[1] = 10
#4. growth rates
R1 = R2 = 2
#5 and 6.  simulate and plot!
for(i in 2:T)
{
	#get d.d. seed values
	seeds = compet(N1[i-1],N2[i-1])
	N1[i] = R1*seeds$N1
	N2[i] = R2*seeds$N2	
}
# N1 wins handily
# does either species drive drive the other into extinction?
plot(N1, type="b")
points(N2, type="b", col="red")
legend(40,60, legend="N1 -> black (K=150)\n N2 -> red (K=100)\n")


### 2 
x = y = seq(1,10)

# set up grid of all potential i-to-j transverses
xy = expand.grid(x,y)

# calculate spacial distances by taking the squared differences in the x-coords
# and adding them to squared dist. in y-coords, then taking the sqrt
#
# first outer() function - gets distances between all xi,xj pairs
# second outer() function - gets distances between all yi,yj pairs
# add them together and take the square root
dmat = sqrt(outer(xy[,1],xy[,1],"-")^2 + outer(xy[,2],xy[,2],"-")^2)


# calculate two redistribution matrices
# one for each species
# sp1 having the 1/9 dispersal (1/9th of all dispersers)
ksp1 = ifelse(dmat <= 1.5, 1/9 , 0)
ksp2 = matrix(1/100, nrow=dim(dmat)[1], ncol=dim(dmat)[2])
diag(ksp1) = 1/9
diag(ksp2) = 1/100

# 9.) matrices to hold simulation data
N1mat = N2mat = matrix(NA, nrow=length(x)*length(y), ncol=T)

# first time step should place 10 "units" of each species at each
# location.
N1mat[,1] = N2mat[,1] = 10

# 10.) spatial analog of pt 5 (above)

for(i in 2:T)
{
	# get number of seeds created at each patch
	# are prev. abundances added to or replaced?
	N1seeds = R1*N1mat[,i-1]
	N2seeds = R2*N2mat[,i-1]
	
	# calc post-dispersal abundances and feed into 
	# compet function
	tmp = compet(N1 = N1seeds%*%ksp1,N2 = N2seeds%*%ksp2)
	N1mat[,i] = tmp$N1
	N2mat[,i] = tmp$N2
	
} # end for loop

# did N1 win again?  
# find the time stamp where there are "no" 
# local occurances of N2 (use 0.00001 as arbitrary cut-off)
min(which(apply(N2mat,2,"mean") <= 0.01))
#it does! so plot

plot(apply(N1mat,2,"sum"), xlab="time", ylab="total abundance in system", type="b")
points(apply(N2mat,2,"sum"), type="b", col="red")
#legend(40,60, legend="N1 -> black (K=150)\n N2 -> red (K=100)\n")


### 3

#new matrices
N1matD = N2matD = matrix(NA, nrow=length(x)*length(y), ncol=T)
N1matD[,1] = N2matD[,1] = 10

#re-use the code
for(i in 2:T)
{
	#get disturbance matrix
	disturbance = rbinom(n=100,1,p=0.40)
	N1seeds = R1*N1matD[,i-1]
	N2seeds = R2*N2matD[,i-1]
	
	#NOTE the change in K1
	tmp = compet(N1 = N1seeds%*%ksp1,N2 = N2seeds%*%ksp2, K1=200)
	#tmp = compet(N1 = N1seeds%*%ksp1,N2 = N2seeds%*%ksp2)
	N1matD[,i] = tmp$N1*(1-disturbance)
	N2matD[,i] = tmp$N2*(1-disturbance)
	
} # end for loop


min(which(apply(N2matD,2,"sum") <= 1))
plot(apply(N1matD,2,"sum"), xlab="time", ylab="total abundance in system", type="b")
points(apply(N2matD,2,"sum"), type="b", col="red")
#legend(40,60, legend="N1 -> black (K=150)\n N2 -> red (K=100)\n")



### last part (reuse a ton of code)

# sequence of p values
pVal = seq( 0.3, 0.9, by=0.01)

# list to hold generation at which 
# N2 or N1 go extinct
Eyear = list( N1 = rep(250,length(pVal)), N2 = rep(250,length(pVal)) )

for( j in 1:length(pVal) )
{
# large number of interations
T = 250
 
#new matrices
N1matD2 = N2matD2 = matrix(NA, nrow=length(x)*length(y), ncol=T)
N1matD2[,1] = N2matD2[,1] = 10

#re-use the code
for(i in 2:T)
{
	#get disturbance matrix
	disturbance = rbinom(n=100,1,p=pVal[j])
	N1seeds = R1*N1matD2[,i-1]
	N2seeds = R2*N2matD2[,i-1]
	
	#NOTE the change in K1
	tmp = compet(N1 = N1seeds%*%ksp1,N2 = N2seeds%*%ksp2, K1=200)
	N1matD2[,i] = tmp$N1*(1-disturbance)
	N2matD2[,i] = tmp$N2*(1-disturbance)
	
} # end for loop
	
	#if sp2 reached 0....
	if( Inf != min(which(apply(N2matD2,2,"sum") <= 10)) )
	{
		Eyear$N2[j] = min(which(apply(N2matD2,2,"sum") <= 10))
	}
	
	if( Inf != min(which(apply(N1matD2,2,"sum") <= 10)) )
	{
		Eyear$N1[j] = min(which(apply(N1matD2,2,"sum") <= 10))
	}
}

#par(mfrow=c(2,1))

plot(pVal, Eyear$N1, ylab="Years to extinction (sp1)",xlab="disturbance values", type="b")
plot(pVal, Eyear$N2,ylab="Years to extinction (sp2)",xlab="disturbance values", type="b")


